Correlations and scaling in one-dimensional heat conduction 
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We examine numerically the full spatio-temporal correlation functions for all hydrodynamic quan- 
tities for the random collision model introduced recently. The autocorrelation function of the heat 
current, through the Kubo formula, gives a thermal conductivity exponent of 1/3 in agreement with 
the analytical prediction and previous numerical work. Remarkably, this result depends crucially on 
the choice of boundary conditions: for periodic boundary conditions (as opposed to open boundary 
conditions with heat baths) the exponent is approximately 1/2. This is expected to be a generic 
feature of systems with singular transport coefficients. All primitive hydrodynamic quantities scale 
with the dynamic critical exponent predicted analytically. 

PACS numbers: 05.10.Ln, 75.40.Mg 



I. INTRODUCTION 

Heat conduction in one dimensional systems is a sim- 
ple example of the general problem of singular trans- 
port coefficients lj. Apart from this theoretical signifi- 
cance, there is increasing experimental relevance due to 
the tremendous advances in nanotube technology 0, 0, 
and the importance of understanding their thermal prop- 
erties. It appears that the thermal conductivity k(L), 
as a function of the system size L, diverges with in- 
creasing L for one-dimensional momentum conserving 
systems 0,|MIH0- RecentlyQ, a hydrodynamic descrip- 
tion has been proposed, which obtained the analytic pre- 
diction of 

k(L) ~ L 1 / 3 (1) 

for such momentum conserving systems 00 ■ Numerical 
results on various systems ja, l| > most recently work on 
Sinai and Chernov's pencase model0] and the random 
collision (RC) modelQ, confirm this prediction. 

Despite the numerical confirmation, there is reason to 
investigate the analytical predictions further. The hy- 
drodynamic description makes detailed predictions about 
the dynamical behavior of all conserved quantities in the 
system: energy, momentum and mass. The heat con- 
ductivity is only a single, and rather indirect, verifica- 
tion of the theoretical picture. In addition, there are 
significant corrections to the scaling form predicted for 
k(L), so that it is appropriate to explore whether they 
are indeed simply corrections to scaling for small L, or 
reflect a departure from the hydrodynamic description. 
Finally, the connection between conductivity and corre- 
lation functions through the Kubo formula [l0|. invoked 
in the analytical treatment , is worth further examina- 
tion. 

In this paper, we study the RC model[9|, which was 
recently introduced for one dimensional transport, and 
numerically obtain the full spatiotemporal form of all 
correlation functions between the conserved hydrody- 
namic variables. We find both propagating and disper- 
sive modes, with a dispersion relation in agreement with 
the analytical result. 



The autocorrelation function of the heat current, which 
can be obtained in terms of the primitive hydrodynamic 
quantities has a spatiotemporal structure that is not 
amenable to any simple scaling description for length and 
time scales much larger than (about thirty times) the in- 
terparticle spacing and collision time. However, as we 
will argue, the existence of larger corrections to scaling 
for the heat current is perhaps not surprising. More- 
over, the continuity relation connects the heat current 
to the primitive hydrodynamic quantities, so these cor- 
rections must vanish on suitably large length and time 
scales. In accordance with this expectation, the auto- 
correlation function of the total heat current (integrated 
over the entire system), Cjj{t), is much better behaved. 
The integral 



Cjj(t') 



(2) 



scales with t with an exponent close to 1/3, although 
depending on the model parameters chosen L can be as 
large as ~ 10000 before the asymptotic scaling behavior 
is seen. 

The autocorrelation function Cjjif) was measured 
with two different boundary conditions: open bound- 
ary conditions with heat baths at the ends, and periodic 
boundary conditions. In both cases, I(t) looks the same 
at short times, and has oscillatory behavior for t ~ L. 
For open boundary conditions, the oscillations damp out 
rapidly, and therefore asymptote on a time scale ~ L, 
as assumed in the analytical treatment. Thus J(f — > oo), 
which is connected to the conductivity through the Kubo 
formula, scales as ~ L 1 / 3 . On the other hand, with pe- 
riodic boundary conditions, I(t) continues to grow after 
the oscillations set in, and /(oo) ~ L 1 / 2 . This depen- 
dence on boundary conditions, although highly unusual, 
is consistent with the analytical picture: the dynamical 
exponent of 3/2 would normally convert a i 1 / 3 scaling to 
i 1/2 , but with open boundary conditions the cutoff to 
the time integral is set by the modes propagating to the 
heat baths at the ends, i.e. at t of O(L). 
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II. EARLIER RESULTS 

A large number of systems with dimension d = 1 have 
been analytically^, 03 E3 or numerically^, IE IE [3 
IT3 | found to have a singular heat conductivity, «(!/) ~ 
L Q with a > 0. Various values of a have been obtained 
for different models. Positive values of a are seen cither 
when a model is integrable, when the value of a depends 
on the details of the model, or when the interparticle 
interactions in the model are momentum conserving. 

For any non-integrable system, it is possible to con- 
struct a hydrodynamic description, where it is assumed 
that local thermal equilibrium is reached^, Since 
long range equilibrium is precluded in d = 1, the hydro- 
dynamic description is only in terms of conserved quanti- 
ties: energy, number and (when momentum is conserved) 
momentum densities. This is the standard hydrodynamic 
treatment of a normal fluid. There are two propagat- 
ing sound modes, and a diffusive mode for heat trans- 
port. The addition of thermal noise to the hydrodynamic 
equations leads to singular corrections when d < 2. The 
spreading of a sound pulse is now superdiffusive: the 
width I is related to time t as t ~ Z 2 , with a dynamic 
exponent z exactly equal to 3/2 for d = 1. The heat 
diffusion mode is no longer diffusive: contributions to 
the heat current that are propagating and bilinear in the 
primitive hydrodynamic densities, that are inconsequen- 
tial for d > 2, dominate (for long wavelength and low 
frequency) for d = 1. The autocorrelation function of the 
heat current decays as Cjj(t) ~ t~ 2 ^ 3 . If one integrates 
this function, with open boundary conditions and heat 
baths, the propagating nature of the modes cuts off the 
integral I(t) at t of O(L). Using the Kubo formula, one 
obtains k(L) ~ L a with a — 1/3. 

III. THE MODEL 

There have been numerical studies of heat conduction 
in various one dimensional models @, IE IE El El- O ne of 
the simplest such models is that of a chain of point par- 
ticles undergoing one dimensional elastic collisions, with 
heat baths at the two ends of the chain. The drawback 
of this model is that it is not chaotic 5]; in fact, when 
the particles in the chain are all of equal mass, the model 
is integrable. As explained in the previous section, the 
hydrodynamic description that yields a universal value of 
a = 1/3 is only valid if the system reaches local thermal 
equilibrium. Even when the masses of the particles in 
the chain alternate, a is found to converge very slowly 
towards its asymptotic value, requiring system sizes of 
~ 16, 000 or larger. 

To alleviate this problem, we recently introduced the 
random collision (RC) model 9]. This can be viewed as 
the limit of rough particles confined to a narrow tube, 
with the width of the tube sufficiently small that par- 
ticles cannot pass each other, in the limit that the size 
of the particles (and the width of the tube) approach 




FIG. 1: Grayscale plot of the autocorrelation function of the 
momentum density, C vv (x,t). The vertical and horizontal di- 
rections are x and t respectively. Periodic boundary condi- 
tions are used for the x direction. The length of the system 
is L = 512. The bright regions correspond to intensity peaks, 
and clearly show the propagating sound modes. 

zero. In this limit, the particles move along a one di- 
mensional line, but possess longitudinal as well as trans- 
verse momentum, both of which contribute to the parti- 
cle energy. Interparticle collisions are elastic. Owing to 
the roughness of the particles, particles emerge from a 
collision with random momenta (respecting energy mo- 
mentum conservation and detailed balance). The trans- 
verse degree of freedom effectively serves as a random 
source/sink of longitudinal kinetic energy in any collision, 
while strictly maintaining conservation of total energy. 
The extra randomization introduced by the transverse 
momentum is expected to equilibrate the system much 
more effectively. Numerically, the conductivity k(L) is 
found to scale much better, with L a scaling observed for 
L ~ 100 to 1000. However, depending on the particle 
masses in the model, there is a slight variation in the 
measured a, with a — 0.28 when all the particles have 
the same mass, and a = 0.33 when the particle masses al- 
ternate with a mass ratio of 2.62. This indicates that cor- 
rections to the asymptotic scaling form, although much 
reduced, have not been completely eliminated. 



IV. SPATIOTEMPORAL CORRELATIONS 

In this section, we present numerical results for the 
spatiotemporal dependence of the correlation functions 
of the primitive hydrodynamic densities: energy, momen- 
tum and number. The motivation for carrying out this 
study has been discussed at the beginning of this paper. 

Figure^shows the autocorrelation function of the mo- 
mentum density, C vv (x, i), as a function of position x and 
time t. The size of the system is L = 512, all particles 
have unit mass, and periodic boundary conditions in x 
are used. One can clearly see two pulses propagating in 
opposite directions at constant speed, corresponding to 
the two sound modes. Even when the pulses collide with 
each other, they emerge essentially unaffected. The fig- 
ure also shows a slight broadening of the pulses as t is in- 
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FIG. 2: Autocorrelation function of the momentum density, 
C vv (x,t), as a function of x for t = 32, 64 and 128. Periodic 
boundary conditions are used for x, and the length of the 
system is L = 512. Since C vv is symmetric in x, only x > 
is shown. As is shown in the figure, the data fits a Gaussian 
form excellently, propagating outwards at constant speed and 
broadening as ~ i 2//3 . 



FIG. 3: Crosscorrelation function of the mass and momentum 
density, C M (i, t), as a function of x for t = 32, 64 and 128. 
Periodic boundary conditions are used for x, and the length of 
the system is L = 512. Only x > is shown, since the function 
is antisymmetric. The large propagating peak, corresponding 
to the sound mode, is fitted to a Gaussian in the figure. The 
bump near the origin comes from the heat mode, and has not 
been fitted. 



creased. To make this picture quantitative, Figure[21plots 
C vv (x,t) as a function of x for various values of t. The 
lines going through the data points are fits to the func- 
tional form C vv (x,t) = (A/t 2 / 3 )exp[~B(x - ci) 2 /* 4/3 ], 
with A, B and c as fitting parameters that are indepen- 
dent of t |16| . The fitting form works surprisingly well. 
Although for spatial dimension d > 2 one indeed ex- 
pects a propagating Gaussian pulse, for d < 2 the scaling 
theory only predicts that the pulse should spread out as 
~ t 2 / 3 , and does not require a Gaussian form. Similar re- 
sults are shown in Figure [3] for the cross correlation func- 
tion of the number and momentum densities, C nv (x,t). 

Figure 0] shows the autocorrelation function for the 
number density C nn (x,t). In addition to the propagat- 
ing parts, one can see a peak near the origin. This peak 
comes from the heat transport mode, which is diffusive 
for d > 2. The same structure is observed in the auto- 
correlation function for the energy density C ee (x,t). In 
order to extract the heat transport mode, suppressing 
the sound modes, in FigureElwe plot C qq (x,t), the au- 
tocorrelation function of 



q(x,t) — e(x,t) — (h/n)n(x,t), 



(3) 



where h and n are the spatially averaged enthalpy and 
number density respectively. At long wavelengths and 
low frequencies, this corresponds to the heat transport 
mode |15). The peak broadens with time with a width 
~ t 2 / 3 , but the fit is not as good as for the previous plots. 

The scaling of C qq (x,t) requires careful consideration 
for d < 2. The hydrodynamic equation for q(x,t) is 

d t q + V • (vdq + vaSe) = k V 2 T + 0(v 2 ) (4) 

where Kq is the bare heat conductivity, T is the temper- 
ature, 8q, 5e are the deviations of q, e from their average 
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FIG. 4: Autocorrelation function of the mass density, 
Cnn(x,t), as a function of x for t = 32, 64 and 128. Peri- 
odic boundary conditions are used for x, and the length of 
the system is L = 512. The function is fitted to the sum of a 
propagating and diffusing Gaussian, with the width of both 
scaling like ~ t 2 ^ 3 .. The second peak does not fit as well as the 
first one; as discussed in the text, the corrections to scaling 
for this peak are expected to be much stronger. 



values, and a — h/e — 1. (For the RC model, a = 1.) 
When fluctuations around average values are small, the 
second term on the left hand side of Eq. Q is second or- 
der in the fluctuating fields and can be neglected, yield- 
ing diffusive transport for q. However, for d < 2, one can 
see Q that this advective contribution to the heat cur- 
rent, vSq + vaSe, dominates at long wavelengths and low 
frequencies. Thus one expects q to propagate, with the 
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FIG. 5: Autocorrelation function for the heat transport mode, 
Periodic boundary conditions are used for x, and the length 
of the system is L = 512. Comparing with C nn (x,i) one can 
see that the sound peak has been eliminated. The fit shown in 
the figure is to a diffusing form: a Gaussian that is centered at 
the origin, spreading out as t 2 / 3 . The fit is only approximate, 
as was the case for the peak at the origin for C„ n inFigure\$\ 




FIG. 6: Grayscale plot of the autocorrelation function of the 
heat current, Cjj(x,t). The vertical and horizontal directions 
are x and t respectively. Periodic boundary conditions are 
used for the x direction. The length of the system is L — 512. 
The intensity peaks move out linearly from x = as a function 
of time. 



propagating behavior becoming increasingly important 
at long length and time scales. This propagating behav- 
ior was used in Ref. Q to argue for a cutoff of t ~ L 
for correlations in a system of size L with heat baths at 
the ends. We will discuss this issue further in the next 
section. 



V. HEAT CURRENT 

We have seen in the previous section that the con- 
served hydrodynamic quantities satisfy the predictions 
of the scaling theory: two propagating modes which dis- 
perse according to Sx ~ i 2 / 3 , and a heat transport mode. 
We have also seen that the third mode, unlike for d > 2, 
is not a simple diffusing mode. We now turn to the au- 
tocorrelation function of the heat current. For momen- 
tum conserving systems, it can be shown on very general 
grounds 01 that this is the correlation function that is 
related to the thermal conductivity through the Kubo 
formula |To| . 

Figure shows a grayscale plot of the autocorrelation 
function of the heat current, Cjj(x,t), similar to that for 
C vv shown in the previous section. The current j(x, t) is 
defined as 



j(x,t) = y^ 5(xj - x)Vi[ei - (h/n)} 



(5) 



where the sum is over particles and ej is the energy of the 
i'th particle. j(x,t) is the current that corresponds to q; 
its hydrodynamic form can be obtained from Eq.Q as 




50 100 150 

x 

FIG. 7: Autocorrelation function of 
Cjj(x,t), as a function of x for t = 32, 
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64 and 128. Peri- 



j(x, t) = v[5q + Se] - k VT + 0{v 2 ). 



(6) 



odic boundary conditions are used for x, and the length of 
the system is L — 512. The peaks decay, broaden, and shift 
outwards as a function of time. No simple scaling form fits 
the data. 



The propagating nature of the heat current can be seen 
by following the intensity peaks in the figure; the range 
in time over which this can be seen is not as large as for 
C vv , because Cjj decays much faster. As in the previ- 
ous section, Figure [7| shows Cjj as a function of x for 
different values of t. Although the peaks shift outwards 
and broaden with time, unfortunately no clear scaling 
of any form can be seen at small t. This result might 
seem disturbing, since the values of t in Figure though 
not large, are sufficient to see good scaling in the basic 
conserved quantities, as shown in the previous section. 
However, as seen in Eq.©, the hydrodynamic form of 
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FIG. 8: Log log plot of I{t)/L x ^ versus t/L for systems of 
different lengths, where I(t) is defined in Eq/(|5J. The particle 
masses alternate between 1 and 2.62. Open boundary condi- 
tions with heat baths at the two ends are used. The data 
show a good scaling collapse, indicating that I(oo) ~ L 1 / 3 . 

the current has advective contributions that arc bilinear 
in the primitive conserved quantities, which dominate at 
long length and time scales, and other dissipative con- 
tributions as well. In a diagrammatic field theoretic ex- 
pansion, one expects relatively larger contributions from 
short wavelength modes, aggravated by the fact that al- 
though the first term on the right hand side of Eq.© 
dominates the asymptotic scaling behavior of j(x,t), it 
is second order in fluctuations and will therefore have 
large corrections to scaling 18]. (The correlation func- 
tion Cjj [x, t) decays so quickly with t that it is impossible 
to measure it for much larger t than shown in the figure.) 

Since we have argued that the lack of scaling seen for 
Cjj(x,t) is a result of strong subleading corrections to 
its asymptotic behavior, it is reasonable to examine its 
spatial integral, i.e. Cjj at zero wavevector. An added 
motivation is that, if J(t) = f dxj(x,t), the Kubo for- 
mula obtains the conductivity as 

K = f - dt'Cjj{t') (7) 
L Jo 

where (3 is the inverse temperature. From the definition 
of Eq.lO, the right hand side of this equation is (3 2 I(t). 
We first show the results for I(t) for somewhat different 
conditions than those used so far in this paper. 

Figure |H1 is a log-log plot of /(^/L 1 / 3 as a function 
of t/L for different values of L. Unlike the data shown 
before this, the system does not have periodic bound- 
ary conditions, but is terminated by heat baths (both 
at unit temperature) at both ends. Also, the masses of 
the particles are not all equal, but alternate between 1 
and 2.62 [l^. The figure shows a power law rise for I(t) 
with exponent i 1 / 3 , which saturates at t ~ L. This is in 
accordance with the analytical prediction Q. There are 
a few damped oscillations at integer multiples of t ~ L, 



FIG. 9: Log log plot of I(t) versus t for L = 16384, where 
I(t) is defined in Eq/(|5J. All particles have unit mass. The 
slope shows an upward drift, from 0.24 (from t = 10 to t = 
100) to 0.29 (from t = 400 to t = 20000), consistent with 
an eventual asymptotic slope of 1/3. The figure also shows a 
similar plot for a chain whose masses alternate between 1 and 
2.62, shifted downwards for clarity. A fit to a slope of 0.33 is 
shown; no upward (or downward) drift in the slope is seen. 
Periodic boundary conditions are used; over the range fitted, 
the results are the same with open boundary conditions. 



presumably corresponding to successive reflections from 
the ends of the system. 

When all the particles in the chain are taken to have 
equal mass, there is a slight discrepancy from the ana- 
lytical prediction. Figure |5] shows a log-log plot of lit) 
for this case. The power law growth of I(t) has an ex- 
ponent close to 0.25, but the saturation at t ~ L is still 
very clear. Taken at face value, this exponent would 
imply a conductivity exponent of a = 0.25, in contra- 
diction to the analytical prediction. However there are 
two reasons why this must be regarded simply as lack of 
complete convergence to the asymptotic form. Firstly, 
as seen in this paper, the dynamics of the primitive hy- 
drodynamic variables do agree with the analytical pre- 
dictions; a fit with Sx ~ £ - 75 does not work. Secondly, 
and more directly, simulations on much larger systems, 
with L = 16000, show an upward drift of the slope of 
ln/(i) versus lni, towards 1/3. No such drift of the slope 
is seen for the unequal mass case shown in Figure[5] Note 
that there are also deviations from the asymptotic value 
of a for the equal mass chain when the conductivity is 
measured directly 0, which was in fact the motivation 
for the detailed spatiotemporal measurements reported 
in this paper. 

The most unusual result for the the heat current au- 
tocorrelation function is obtained for L{t) for the case of 
periodic boundary conditions. Figure ITUI is a log- log plot 
of 7(<)/L 1 ^ 3 as a function of t/L. As was the case with 
heat baths at the ends, L(t) grows initially with the form 
~ i 1 / 3 , has damped oscillations at t ~ L, and saturates 
for large t. However, unlike the case with heat baths, 
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FIG. 10: Log log plot of I(t)/L 1/3 versus t/L for systems of 
different lengths, where I(t) is defined in Eq/JS). The parti- 
cle masses alternate between 1 and 2.62. Periodic boundary 
conditions are used. For small t, I(t) is the same as for open 
boundary conditions (shown in Figure |5J. However, beyond 
the oscillatory regime at t ~ O(L), I(t) keeps rising, and I(oo) 
does not scale as ~ L 1 ^ 3 . 
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and the eventual asymptotic value scales approximately 
like L 1 / 2 . In hindsight, this result is reasonable: with- 
out the heat baths, the dynamics are not cut off at the 
ends of the system, and the dynamic scaling exponent of 
z = 3/2 saturates the t 1 ' 3 growth at ~ L 1 ! 2 . Alterna- 
tively, the saturation of I(t) is set by the time a pulse 
takes to spread across the system rather than propagate 
across it. It is not clear exactly why the spreading cuts 
off the growth of I(t), though from a scaling viewpoint 
one could argue that in the absence of the propagation 
timescale, the only possibility is L ~ i 2 / 3 . 

VI. DISCUSSION 



We have seen that the dynamics of the primitive hydro- 
dynamic quantities, mass, momentum and energy den- 
sity, fit a scaling description with the predicted dynamic 
exponent. The scaling of the heat current autocorrela- 
tion function is more problematic. However, we have 
argued that this is to be expected, because the asymp- 
totically dominant part of the current operator is bilinear 
in the primitive quantities and has a small bare value in 
the sense of the renormalization group compared to cor- 
rection terms. The same slow convergence to asymptotic 
scaling that afflicts the heat current autocorrelation func- 
tion has been noticed for the thermal conductivity and 
remarked on by various authors. However, the results 
for the primitive hydrodynamic quantities make it clear 
that this is indeed simply a case of slow convergence to 
the asymptotic limit rather than a discrepancy with the 
analytical prediction. Most remarkably, the large time 
limit of the spatiotcmporal integral I{t) of the heat cur- 
rent autocorrelation function scales with L as ~ L 1 / 3 for 
a system terminated with heat baths, and as ~ L 1 ! 2 for 
periodic boundary conditions. This implies that in ap- 
plying the Kubo formula to systems with singular trans- 
port coefficients one must be careful about the boundary 
conditions being used. 



FIG. 11: Log log plot of /(oo) versus L for open and periodic 
boundary conditions. The straight lines show fits to ~ L 1 ^ 3 
and ~ L 1 / 2 respectively. 
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